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Abstract 

One of the most efficient ways to produce unconditional simulations is with the spectral method using fast Fourier transform 
(FFT) Q. But this approach is not applicable to arbitrary surfaces because no regular grid exists. However, points on the arbitrary 
surface can be generated randomly using uniform distribution to replace a regular grid. This paper will describe a nonstationary 
kernel convolution approach for data on arbitrary surfaces. 

Index Terms 

Kernel Convolution; Turning Band; Nonstationary Simulations; Large Data 

I. Introduction 

The approach to producing unconditional simulations for arbitrary surfaces (sphere, ellipsoid, geoid, torus, etc.) will be 
described in this paper. The approach is also applicable to nonsmooth surfaces. This paper is based on ideas from a turning 
band approach (Q, Q, 0, and 0) and kernel convolution (0, |[7|, and 0). Similar to the turning band, all processing is 
done on a set of random lines. The difference is in using a random set of points covering an area of interest and projecting 
them onto the set of random lines. For example, the set of uniformly distributed points on the surface of a sphere is projected 
onto the set of random lines. The kernel convolution is applied to each line, and the reconstruction of the process is performed 
similarly to the turning band approach. Using different kernels will produce a nonstationary solution. 

II. Kernel Convolution by Turning Band 

Here are the steps of the algorithm: 

1) From a uniform distribution, simulate N points on a surface. 

2) From a standard normal distribution, simulate value for each point. 

3) Simulate random or pseudorandom directions (the algorithm to simulate pseudorandom directions can be found in 0 
and 0). 

4) Cycle over all directions. 

a) Project points to a direction. 

b) Apply kernel convolution (this is usually done using FFT). 

5) For each location where the result is desired, project that location to all directions and average the results of kernel 
convolutions. 

The applied one-dimensional kernel kl{^) will correspond to a three-dimensional kernel ks{^) through integral 

^3 ^ fci ^ dn, (1) 

where (x, y) is a scalar product. 

The kernel distance is the direct distance in three-dimensional space. 

This is similar to the turning band method, and the relationship between three-dimensional or two-dimensional covariance 
with one-dimensional covariance can be found in 0, 0, and 0. The same relationship holds true for the kernel functions. 
In a three-dimensional case, the formula is 

ki{h) = {h^ks{h)y. (2) 

^Replacing uniform distribution on the surface by pseudorandom uniform distribution will produce an even coverage and guarantee the stability of the 
algorithm. 


In two dimensions for kernel k 2 {^) 


tl ^ 2 

ki(h) = k2{^) ^ h f — ^ dx = ^2(0) + h f sin ( 0 )) d 0 . 


0 ^ ^ 0 

III. Kernel Convolution for Unconditional Simulation 


Applying the algorithm from the section ‘ Kernel Convolution by Turning Band ’ for any location s will produce a random 
field approximately equal to 


N 


(k3(Si - ^ ■ ^i) . 


(3) 


i=0 


The approximation is due to the finite number of directions. 
The variance at location s is approximately equal to 


N 


Ym-s). 


i=0 


(4) 


For the sphere, the estimate of 0 does not depend on location s. This does not hold true for other surfaces. However, 
for any surface with finite sampling N, it is not a constant. Calculating this value is important in standardizing unconditional 
simulations. The estimate of can be found by applying the algorithm from the section ‘ |Kernel Convolution by Turning] 
Band ’ for the square of the kernel for the same set of points with values of one. Dividing the approximation of Q by the 
square root of the approximation of Q will produce an approximation of the Gaussian process with unit variance. 

Another way to estimate Q is by producing many unconditional simulations and estimating the variance from them. This 
will also allow the estimation of covariances. 

Compared to the turning band approach, using this algorithm is significantly slower. The slowest part is projecting simulated 
points to directions. Notice that the simulated points do not depend on the kernel. Therefore, to speed up the approach, it is 
possible to simulate points and project them to directions in advance. Combining them by bins will significantly reduce the 
amount of memory necessary to store precalculated results. 


IV. Conditional Simulations 

To perform conditional simulations, it is necessary to know the covariance function. For a sphere, the covariance function is 
isotropic and does not depend on location. Therefore, it can be derived directly from the kernel as described in 0 and |TQ| . 
For arbitrary surfaces, it is possible to make many unconditional simulations and from them estimate the covariance function. 
Conditional simulations are obtainable from unconditional simulations through solutions of the Kriging system (TT), (^, and 
(T3) (13 describes the approach to construct conditional simulations without abrupt changes. 

V. Class of Exponential Kernels 

Let the kernel function have the form 


where Pi is a polynomial. 

The result of the convolution can be obtained by the sequential application of the kernel. The technique is described in (H) 
and |T^ with the difference being that the integration is replaced by summation. This reduces the complexity of the kernel 
convolution from 0{n • Inn) to 0{n), where n is the number of bins on the line. 

For finite kernel 


ksih) = < 


Y 



if h<T, 


( 6 ) 


0 


otherwise, 


where T is the kernel width. To be able to evaluate ki{h) from 0, it is necessary to satisfy lim ks{h) = 0. Otherwise, 
ki{h) will not be finite. 

There are several advantages of using the finite kernel: 

• Ability to produce unconditional simulations for a large area by partitioning it into small pieces (tiles) for independent 
processing. 

• Ability to have a fiexible class of kernels in form ([^. 
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The ability to have a flexible class of kernels is important for nonstationary simulations. It allows not only changes in kernel 
bandwidth but in their shape as well. One way to construct a flexible class of kernels in form ([^ is by changing the polynomial 
part. It is reasonable to use decreasing kernels because the majority of physical processes gradually decreases dependence over 
distance. From (13 a Bernstein density (a weighted sum of beta probability distribution functions defined on the unit interval 
[0,1]) is defined as 

m m 

p{x \ m^w) = ^ (^Wj • Beta{x | j,m — j + l)^ , > 0 , ^ Wj = 1 ^ 

3=1 ^ 3=1 


where m is the number of distribution functions, w are weights, Beta{^) is the beta probability density function 

ml 


Beta{x I m — j + l) = — 


Forming 


/ 1 


• ^ • (1 — x)^ ^ 


/ III' p 

p(^x \ m^w) = Wj / Beta{y \ j^m — j 1 ) 

7 = 1 ^ 

X \ X 


dy 


and multiplying by produces a flexible family of kernels. The choice of 0 is arbitrary. If the value is too large, the floating 
point evaluation will lead to a high round off error. If the value is too small, it will not produce good family of kernels. Let’s 
define 0 = 3. 

In the case of m = 2, this becomes 


Wl (1 — d- W2 ~ 


_ h 

e 3 


0, 


if < 1, 
otherwise, 


where 0<u;^<l,i = l,2, and Wi d- W 2 = 1. They are shown in Figure [^. 


(7) 




a) b) 


Fig. 1. Example of kernels 0 mi 


i 

15’ 


m 2 = 1 — mi, ^ = 0..15. a) With exponent, b) Without exponent. 


In the case of m = 3 

(1 — h)^ + UJ 2 (1 — h)‘^ (1 + 2h) d-w^i^l — 

0 , 

\ 

where 0 < < 1, i = 1, 2, 3, and Wi d- W2 d- = 1. They are shown in Figure]^. 


if h < 1, 
otherwise, 


( 8 ) 
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From ([^, it follows that ki{h) has the same form as ([^. Notice that kl{h) also has the same form, and the one-dimensional 
kernel corresponding to kl(h) will be of the same form. 

It is possible to extend this approach for anisotropic kernels by adjusting weights and one-dimensional kernel shape in each 
direction or by stretching the space. 


VI. Examples 


Figure shows unconditional simulation on a torus. 



Fig. 3. Unconditional simulation on a torus with inner radius 1, outer radius 3, and kernel e Number of simulated points N = 67,108, 864. Number 
of random directions = 1, 024. 
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Figure 


shows nonstationary unconditional simulation on a spheroid with eccentricity equals -. The kernel changes smoothly 


in shape Tfom left to right. 1,024 unconditional simulations were used to estimate variance and standardize unconditional 
simulations. Due to changes in the kernel, the right and left parts of the spheroid have different covariance structures. This is 
visible as a smoother random field on the right compared to the random field on the left. 



Fig. 4. Unconditional simulation on an oblate spheroid with semimajor axis equals 2, and semiminor axis equals 1. The surface is shown as elevation to 
depict changes in smoothness. The kernel has the form of (^. It changes smoothly from the kernel on the left • e~^ to the kernel on the right 

(1 — 3h)^ ■ e~^ with bandwidth equals -. Number of simulated points N = 262,144. Number of random directions = 1, 024. 


VII. Convergence Issue 

When the kernel bandwidth tends to be small, the convergence of ^ to the true kernel O becomes slow. The number of 
required directions makes this method inefficient. This is because the distribution of the number of points that are farther away 
from the center of the kernel (location s in ([^) in three dimensions increases squarely. This is even true for some directions 
when points are only located on the surface. 


VIII. Tile Approach 

Finite kernel and overlapping tiles overcome the issue described in the previous section. For example, in two dimensions, if 
all tiles are 2 by 2 and start from integer numbers (corners can be rounded), then the kernel with a bandwidth not exceeding 

i will always be completely inside at least one tile; see Figure This is applicable in any number of dimensions. 

The bandwidth should not be significantly smaller than the tile width, and the exponential part in the kernel can be dropped. 
The fiexible family of kernels without the exponential part is shown in Figures and |^. 
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Fig. 5. The red circle is the area covered by a finite kernel. The blue rounded square is the tile completely covering the finite kernel. 


Figure shows the effect of using the tile approach. The small discontinuities are visible in the horizontal and vertical lines 
passing through the center. Increasing the number of bands will reduce this effect. Figure shows the result of evaluating 
^ and 0 directly. 



Fig. 6. Unconditional simulation in two dimensions by kernel k 3 {h) = (1 — 2h)^ with bandwidth equals -. Number of simulated points N = 65, 536. 
a) The tile approach with the number of random directions = 1, 024 (uniformly distributed in three dimensions), b) Direct calculation. 


If a continuous result is required, overlapped tiles can be joined smoothly. This will increase calculations by the number of 
tiles necessary to process. One way to reduce the number of tiles is to use a triangular grid (use nonequilateral triangles for 
arbitrary surfaces). For each node of the triangular grid, define a tile that covers all neighboring triangles and at least a kernel 
bandwidth. For each triangle, the barycentric coordinate system will give weights for the tiles. There will only be a maximum 
of three tiles with nonzero weights; therefore, the complexity is increased three times. 

IX. Using Integer Arithmetic 

In the turning band approach, the directions are chosen randomly or pseudorandomly. In this section, another approach will 
be described to make directions proportional to integer numbers. Using such directions will remove errors from applying fixed 
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steps (bins) in each direction. 

Let’s take a sphere of radius R = VS, where S' is a natural number. If the radius R approaches infinity, all integer vectors 
inside the sphere (all integer vertices except the center of a coordinate system) tend to be uniformly distributed. However, the 
unique directions are not uniformly distributed (see example in Figure |7^). But assigning a proper weight to each direction 
will resolve it. The appropriate weights can be assigned as proportionate to the number of vectors with the same direction or 
to the area of the spherical Voronoi diagram cells. 



a) b) 

Fig. 7. a) All integer directions are shown as dots on the sphere for S = 100. The number of unique directions equals 1, 729 (counting v and —v once), 
b) A uniform set of integer directions satisfying that any two directions are at least 3° apart. The number of directions equals 1, 405. 



Selecting only directions separated by some degree will make a uniform coverage (see Figure |^). The list of integer 
directions is shown in the appendix This is similar to using a set of pseudorandom directions as described in |[8| and Q. 

Because each direction is represented as a - {x^y^z), where ct is a real number (square root of a natural number) and x, y, 
and z are integers, any coordinate with integer components will be projected to a position a • i, where i is an integer number. 

If random points are placed on an integer grid, there are no errors due to binning. If the resultant location also has an integer 
coordinate, there is no additional error. Notice that the kernel 0 or 0 can be calculated for any position on direction, see 
(Tg and p^ . The turning band approach 0, g, g, and 0 might also benefit from using integer arithmetic. 

X. Conclusion 


The new approach to producing unconditional nonstationary simulations on arbitrary surfaces is equivalent to the kernel 
convolution approach with the kernel shape depending on location (|[^, Q, and Q)- Applying a kernel to the arbitrary surfaces 
produces a larger class of valid covariance functions than the class of valid covariance functions in three-dimensional space. 

Precalculating projections of random points in all directions for all tiles might be necessary for the efficiency of the approach. 
Nevertheless, even with precalculations, this approach is not efficient. The better approach is based on kernel convolution using 
FFT as described in (T^. 

New flexible classes of kernels constructed from polynomials are described in this paper. Because they are linear combinations 
of some basis kernels, constructing covariances for all pairs of basis kernels is sufficient to calculate covariance between any 
two kernels. Unconditional simulation for a nonstationary kernel is a weighted sum of unconditional simulations corresponding 
to basis kernels for the same set of random points. 

Using integer arithmetic, an error introduced in the turning band approach due to binning along each direction can be 
avoided. 
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Appendix. List of Integer Directions 


Integer 

directions that are 

at least 2 

1° apart: 

(1,0,0), 

(1,1,0) 
A ’ 

(1,1,1) 
A ’ 

(2,1,0) 
A ’ 

(2,1,1) 
A ’ 

(3,1,0) 

Ao ’ 

(2,2,1) 
3 ’ 

(3,1,1) 

AT 

(3,2,0) 

(4,1,0) 

(3,2,1) 

(4,1,1) 

(3,2,2) 

(3,3,1) 

(4,2,1) 

(4,3,0) 

(5,1,1) 

(5,2,0) 

(6,1,0) 

(3,3,2) 

(4,3,1) 

A3 ’ 

VT7 ’ 

Ai ’ 

As ’ 

A7 ' 

’ A9 

’ AT 

’ 5 ’ 

A7 ' 

’ A9 

’ A7 

’ A2 ' 

’ A6 

(5,2,1) 

(4,3,2) 

(4,4,1) 

(5,2,2) 

(5,3,1) 

(6,2,1) 

(7,1,1) 

(4,3,3) 

(5,3,2) 

(5,4,1) 

(6,3,1) 

(9,1,0) 

(5,4,2) 

Ao ’ 

Ao ’ 

A3 ’ 

A3 ’ 

As ^ 

’ AT 

’ AT 

’ AT ’ 

As 

’ A2 

’ A6 

’ A2 ^ 

’ As 

(6,3,2) 

(6,4,1) 

(7,3,1) 

(7,4,0) 

(8,2,1) 

(5,4,3) 

(7,3,2) 

(10,1,1) 

(6,4,3) 

(6,5,2) 

(6,6,1) 

(7,4,2) 

(7,6,0) 

7 ’ 

’ 

a/m ’ 

As ’ 

Ao ’ 

^Ap ’ 

A2 ’ 

Ap2 

’ AT 

’ As 

; A3 

’ A9^ 

’ A5 

(10,2,1) 

(5,5,4) 

(7,6,1) 

(10,3,: 

1) (7,4 

,4) (8,€ 

1,1) (9,: 

5,1) (11, 

,3,2) (11,4,1) 

(8,7,2) 

(11,5,1) 

(7,6,5) 

\/T05 

As 

’ A6 

’ ATo 

: ’ 9 

’ Aoi ’ A07 ’ A34 ’ 

A38 ’ 

ATT ’ 

ATT ’ 

ATo 

(11,4,3) 

(7,6,6) 

(8,7,4) 

(10,5,4) 

(11,7, 

1) (17,1 

,1) (18, 

1,0) (13, 

4,3) (16,4,1) (18,3,1) ( 

10,7,6) 

(10,9,4) 


■JtM ’ 11 ’ ^/V^ ’ JlW ’ ^/Y^\ ’ VM. ’ 

(11,11,1) (12,8,3) (12,7,6) (14,12,1) (17,12,1) 

’ Vm ’ ’ Vui ’ ViM 

each vector all signs and permutations of its components. 


’ VlM ’ Vf73 ’ \/3M ’ ’ VWf ’ 

The complete list of directions is formed by considering for 
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